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Abstract 

Dynamics of an ensemble of TV-unit FitzHugh-Nagumo (FN) neurons sub- 
ject to white noises has been studied by using a semi- analytical dynamical 
mean-field (DMF) theory in which the original 2iV-dimensional stochastic 
differential equations are replaced by 8-dimensional deterministic differential 
equations expressed in terms of moments of local and global variables. Our 
DMF theory which assumes weak noises and the Gaussian distribution of 
state variables, goes beyond weak couplings among constituent neurons. By 
using the expression for the firing probability due to an applied single spike, we 
have discussed effects of noises, synaptic couplings and the size of the ensemble 
on the spike timing precision, which is shown to be improved by increasing 
the size of the neuron ensemble, even when there are no couplings among 
neurons. When the coupling is introduced, neurons in ensembles respond to 
an input spike with a partial synchronization. DMF theory is extended to 
a large cluster which can be divided into multiple sub-clusters according to 
their functions. A model calculation has shown that when the noise inten- 
sity is moderate, the spike propagation with a fairly precise timing is possible 
among noisy sub-clusters with feed-forward couplings, as in the synfire chain. 
Results calculated by our DMF theory are nicely compared to those obtained 
by direct simulations. A comparison of DMF theory with the conventional 
moment method is also discussed. 



PACS No. 87.10.+e 84.35.+i 05.45.-a 07.05.Mh 



Typeset using REVTgX 



E-print: cond-mat / 0206135 



t E-mail : hasegawa® u-gakugei .ac.jp 



1 



I. INTRODUCTION 



It has been controversial how neurons communicate information by firings or spikes |]J- 
0. Much of debates on the nature of the neural code has been mainly focused on the 
two issues. The first issue is whether information is encoded in the average firing rate of 
neurons (rate code) or in the precise firing times (temporal code). Adrian first noted the 
relationship between the neural firing rate and the stimulus intensity, which forms the basis 
of the rate code. Actually firing activities of motor and sensory neurons are reported to 
vary in response to applied stimuli. In recent years, however, an alternative temporal code 
has been proposed in which detailed spike timings are assumed to play an important role in 
information transmission: information is encoded in interspike intervals or in relative timings 
between firing times of spikes M- [IU|. Indeed, experimental evidences have accumulated 



in the last several years, indicating a use of the temporal coding in neural systems fll 
|15fl . Human visual systems, for example, have shown to classify patterns within 250 ms 
despite the fact that at least ten synaptic stages are involved from retina to the temporal 
brain ||15|| . The transmission times between two successive stages of synaptic transmission 
are suggested to be no more than 10 ms on the average. This period is too short to allow 
rates to be determined accurately. 

The second issue is whether information is encoded in the activity of a single (or very few) 
neuron or that of a large number of neurons (population or ensemble code). The population 
rate-code model assumes that information is coded in the relative firing rates of ensemble 
neurons, and has been adopted in the most of the theoretical analysis fI6 |. On the contrary, 
in the population temporal-code model, it is assumed that relative timings between spikes 
in ensemble neurons may be used as an encoding mechanism for perceptional processing 



17| - [19| . A number of experimental data supporting this code have been reported in recent 
years [^TJ. For example, data has demonstrated that temporally coordinated spikes can 
systematically signal sensory object feature, even in the absence of changes in firing rate of 
the spikes f22|j . 

It is well known that neurons in brains are subject to various kinds of noises, which 
can alter the response of neurons in various ways. Although firings of a single neocortical 



neuron in vitro are precise and reliable, those in vivo are quite unreliable [23]. This is due to 



noisy environment in viro, which makes the reliability of neurons firings worse. The strong 
criticism against the temporal code is that spikes are vulnerable to noise while the rate code 
performs robustly in the presence of noise but with limited information capacity. It has been 
shown, however, that the response of neurons is improved by background noises against our 
conventional wisdom. The typical example is the stochastic resonance (SR), in which weak 
noises enhance the transmission of subthreshold signals (for review see Refs. [E3 [23])- It 
has been shown that noise of appropriate magnitude linearlizes the response of neurons, 
which leads to SR and maximizes input-output correlation, transformation and coherence 
(for review see Ref. |26j). Recently, it has been demonstrated that noises can enhance the 



firing-time reliability of neurons stimulated by weak periodic and aperiodic inputs 
We may expect that a population of neuron ensembles plays important roles in the response 
of neurons subject to noises. Actually SR in HH neuron ensembles has first demonstrated 
for a single input by Pei, Wilkens and Moss POfl . Subsequently this pooling effect has been 
supported for aperiodic |H]] ]32| and periodic (analog) signals |33|] and for spike-train inputs 
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S3 1 ||35|| . Quite recently, SR for a transient spike signal in large-scale HH neuron ensembles 
has been studied by using the wavelet analysis [Q . It may be possible that the firing-time 
precision is also improved by increasing the size of neuron ensembles. 

A small patch of cortex may contain thousands of similar neurons, each connecting with 
hundreds or thousands of other neurons in that same patch or in other patches. The under- 
lying dynamics of individual neurons includes a variety of voltage dependent ionic channels 
which can be described by Hodgkin-Huxley-type differential equations. Computational neu- 
roscientists have so far tried to gain understanding of the properties of neuron ensembles 
with the use of two approaches: direct simulations and mean-field (MF) theories. Simula- 
tions have been made for large-scale networks mostly consisting of integrate-and-fire (IF) 
neurons. Since the time to simulate networks by conventional methods grows as iV 2 with N, 
the size of the network |36| , it is rather difficult to simulate realistic neuron clusters, in spite 
of recent computer development. In MF theories |37|]- ||40|| , dynamics of globally coupled 



large-scale networks is described by a flow of phases or the population activity, which deter- 
mines the fraction of the firing rate of neurons. The stability condition for synchronous and 
asynchronous states of neuron clusters has been investigated. Quite recently, the population 
density method has been developed as a tool modeling for large-scale neuronal clusters ||41|| - 
pE2|] . In these MF approaches the macroscopic variable of interest is the firing rate, following 
the rate-code hypothesis. However, only little MF approaches have so far proposed based on 
the temporal-code hypothesis |39| . 

The purpose of the present study is to construct a dynamical mean-field (DMF) theory 
based on the temporal code hypothesis, generalizing the method previously proposed by 
Rodriguez and Tuckwell (RT) |[43|| - ||46|| . In the RT theory, the dynamics of the membrane 
potential of a neuron subject to white noises is studied by replacing stochastic differential 
equations (DEs) by deterministic DEs described by moments of state variables. RT's general 
theory has first applied to a single FitzHugh-Nagumo (FN) neuron [43| jPj] and then a 
Hodgkin-Huxley (HH) neuron [3B]. In the case of a single FN neuron, for example, two 



stochastic DEs are replaced by five deterministic DEs, for which an improvement to the RT 
theory has been recently proposed [47]. When the RT theory is applied to a iV-unit FN 
neuron network, 2N stochastic DEs are replaced by N eq = N (2N + 3) deterministic DEs 
f4~3| |46|] . For example, in the case of N = 100, we get N eq = 20300, which is too large 
to perform calculations for neuron ensembles. In their subsequent paper |45[], the result 
of ensemble neurons is transplanted to the Fokker-Planck (FP) equation for the transition 
probability density, which is a partial differential equation with 2N+1 independent variables. 
Solving such a FP equation is a hard computational task. We will present in this paper, an 
alternative MF theory for N FN neuron ensembles, replacing original 2N stochastic DEs by 
eight deterministic DEs which are expressed in terms of means, variances and covariances 
of local and global state variables. 

There are several nonlinear models which have been used for a study of neuron activities. 
Among them we employ here the FN model ESI |49| because it is relatively simple and 



amenable to analysis although the FN model does not have as firm an empirical basis as 
conductance-based model like the HH model. The property of the FN model has been 
intensively investigated. In recent years, SR of a single FN neuron p0| - |52| and FN neuron 
ensembles [[H|] [[32] |53| |]54j have been studied. 

The paper is organized as follows: In Sec. II, we have developed a DMF theory for N 
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FN neuron ensembles, expanding the original stochastic DEs in terms of deviations from 
means to get variances and covariances of local and global variables. We compare our 
DMF theory with conventional RT's moment method ||43|| , showing that the former may be 
derived from the latter. Some calculated results are reported of the response of ensemble 
neurons to a single spike with white noises. It will be shown that the spike firing precision 
is improved by increasing the ensemble size and the synaptic couplings, as expected. In Sec. 
III. DMF theory is extended to a large cluster consisting of multiple sub-clusters and model 
calculations are reported. The final Sec. IV is devoted to conclusions and discussions. 



II. NEURON ENSEMBLES 

A. DMF approximation 

We assume a neuron ensemble consisting of A-unit FN neurons. Dynamics of a single 
FN neuron i in a given ensemble is described by the nonlinear DEs given by 

^ = F[ Xi (t)} - c yi {t) + 4 c) (t) + &\t) + CM, (1) 

= bxi(t) - fa® + e, (i = l-N) (2) 

where F[x(t)} = k x(t) [x(t) - a] [1 - x{t)}, k = 0.5, a = 0.1, b = 0.015, c = 1.0, d = 0.003 
and e = |43| jPJ, Xi and yi denote the fast (voltage) variable and slow (recovery) variable, 
respectively, and is the independent Gaussian white noise with < £j(t) >= and 
< &(*) &(*0 >= Pi $ij s (t - ?), the bracket < ■ > denoting the average over stochastic 
random variables |)5]. In Eq. (1), l[ (t) denotes the coupling term given by 



Ii C \t) = £ E G( Xj (t)), (3) 

where w stands for the coupling strength and G(x) = 1/[1 + exp[— (x — 9)/ai\] is the sigmoid 
function with the threshold 9 and the width a. The self-coupling terms are excluded in Eq. 
(3), where we have adopted the normalization factor to be N~ x instead of (N — 1) _1 for a 
later study of the N = 1 limit. J (e) (t) expresses an external, single spike input applied to 
all neurons, given by 

I^(t) = A Q(t - t^) e(t in + T w - t), (4) 

where Q(x) = 1 for x > and otherwise, A stands for the magnitude, ti n the input time 
and T w the width. 

After RT |4~3| [4~3], we will express these nonlinear DEs by moments of variables. First 



we define the global variables for the ensemble by 

X(t) = (l/N) £ a*®, (5) 

i 

Y(t) = (l/N) E Vi (t), (6) 
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and their averages by 

^(t) = <X(t)>, (7) 

fi 2 (t) = < Y(t) > . (8) 

Next we express the differential equations given by Eqs. (1) and (2) in terms of the deviations 
from their averages defined by 

5xi(t) = Xi (t) -/ii(t), (9) 

6y i (t)=y i (t)- t i2(t), (10) 



to get (the argument t is hereafter neglected) 



and those between global variables, given by 



(11) 
(12) 

(13) 

(14) 
(15) 
(16) 



-cpL 2 - c5 yi + 4 C) + 4 e) + 

^1 — b/i! - dfx 2 + b5xi - dSyi + e, 
at 



with 



4 C) = - ^)G(fii) 

+^ E [G'MSxj + \g^\^)8x) + l -G^{^)8x*\). 

We define the variances and covariances between local variables, given by 

71.1 = ^ E < 6x2 i >. 

i 

72.2 = E < S Vi >' 

i 

7l,2 = ^E <fa > ^ >' 



p M = < 5X 2 >, (17) 

P2,2 = < ^ >, (18) 

Pi, 2 = < 5y >, (19) 

where <5X = — and 5Y = Y(t) — ^(t)- It is noted that 7 Kj a expresses the 

spatial average of fluctuations in local variables of X{ and yi while p K \ denotes fluctuations 
in global variables of X and Y . We assume that (1) the noise intensity j3 is weak and (2) the 
distribution of state variables take Gaussian form. The first assumption allows us to expand 
the quantities in a power series of fluctuation moments around means. As for the second 
assumption, numerical simulations have shown that for weak noises, the distribution of x(t) 
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of the membrane potential of a single FN neuron nearly obeys the Gaussian distribution, 
although for strong noises, the distribution of x(t) deviates from the Gaussian, taking a 
bimodal form (see Fig. 8 of Ref. [33] and Fig. 3 of Ref. [0). Similar behavior of the 



membrane-potential distribution has been reported in a HH neuron model [[33j p8[ . When 
adopting the Gaussian assumption, we may express the average of fluctuations in terms of 
the first and second moments only. It is noted that we impose no conditions on the coupling 
strength. After some manipulations, we get DEs for means, variances and covariances, given 
by (for details see Appendix A): 



dt 
dp 2 
~dt 

dt 

d~?2,2 

dt 

dt 
dpi,i 
dt 

dt 
dpi,2 
dt 



fo + hli,i -cp 2 + w (1 
b/ii - dp 2 + e, 



1 



U + I (e) (t), 



2[(/i + 3/ 3 7i,i)7i,i - c 7 i, 2 ] + 2w(p 1)1 - ^) U x + P 
2(671,2 - ^72,2), 

+ (fl + 3/371,1 - ^)7l,2 - C7 2 ,2 + W (pi, 2 



7l,2 - 



Uu 



with 



2[(/i + 3/371,^1,1 - cp lj2 ] +2w(l- ^) p hl V x + £ 

2(6/71,2 - dp2,2), 

6pi,i + (fi + 3/371,1 - d)p h2 - cp 2 ,2 +w(l-—) Pl, 2 Ux, 



U = g + 92ji,i, 
Ui = gi + 3flf 3 7i,i, 
f e = (l/£\)F^(p 1 ), 
g e = (l/£\)G^(p l ). 



(20) 
(21) 
(22) 
(23) 
(24) 
(25) 
(26) 
(27) 



(28) 
(29) 
(30) 
(31) 



where (3 2 = (1/JV) 

The original 2A^-dimensional stochastic differential equations given by Eqs. (1) and (2) 
are transformed to eight- dimensional deterministic differential equations, which show much 
variety depending on model parameters such as the strength of white noise (/?), couplings 
(w) and the size of cluster (N). 



B. Derivation of DMF theory from RT's moment method 

Before discussing the property of our DMF theory, we will show that it can be derived 



from RT's moment method [43]. In the case of a single FN neuron (N = 1), DMF theory 
agrees with the RT theory as shown in appendix B, where some limiting cases of Eqs. (20)- 
(27) are examined. In the case of FN neuron ensemble (N > 2), the RT theory defines 
means of variables of Xj and yi for the neuron i given by 
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m[ = < Xi >, (32) 
ml = < Vi >, (33) 

and calculate variances and covariances between local variables as given by 

Ci'j = < Ax t A Xj >, (34) 
C% = <A yi A yj >, (35) 
Ci% = < Ax t Ay, >, (36) 

where Axi = Xi — m\ and Ay,i — y% — m\. Variances are given by setting i — j in Eqs. 
(34)-(35). Adopting the same assumptions as our DMF theory: (1) weak noises and (2) the 
Gaussian distribution for state variables, we get DEs for these moments given by (for details 
see Appendix A) 

^ = fo + m:, - cmi + % £ [ 9 i + gidi] + , (3?) 



dm\ 
~df 
dCij 
dt 



bm\ - dm\ + e, (38) 

(fi + men +fi+ mv^r. - c(c$ + c%) + n i3 

w 



+^[ E (9i + 3<&t?)C?;{ + J2(9i+ Sg k 3 C k ^)Ci% (39) 



dC^ 
dt 

2 



dC l { 3 



dt 



*(#*) K+3) 

b(Ci$ + C%) - 2dC%, (40) 
bCH + (ft + WZ - d)CH - cC% 



+^ E (9i + ZdlC^Cil ih J = 1 to N) (41) 



where 



/j = (l/€!)F»(mi), (42) 

^ = (l/€!)(?W(mi). (43) 

Now we derive DMF theory from RT's moment method. We adopt the approximation 
given by 

ml = Aii, (44) 

which yields = and g\ = gi [Eqs. (28), (29), (42) and (43)], and the approximation 
given by 

/ 3 Ci;i ~ /:;-,., 

g 3 Ctf ~ ^71,1- ( 45 ) 
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in Eqs. (39) and (40). We realize that quantities of p K , 7 Kj a and p Kt \ [Eqs. (5)-(8) and 

(14)-(19)] adopted in DMF theory are expressed in terms of m % K and [Eqs. (32)-(36)] in 
the RT theory as follows: 

= /y E m l ( 46 ) 

i 

7„a=4E^ (47) 

i 

^ = ]sEE^ («,A = 1,2) (48) 

Then, we may obtain, from Eqs. (37)- (41), (44)-(48), alternative DEs for p K , 7 K]A and p K _ A 
which are again given by Eqs. (20)-(27). This implies that DMF theory may be derived 
from RT's moment method if we adopt the assumptions given by Eq. (44) and (45) |56| . 

Taking into the symmetry relations: C l ^\ = Cjj\, we get the number of DEs to be 
N eq = 2N + N(2N + 1) = N(2N + 3) in the RT theory [Eqs. (37)-(41)] while N eq = 8 
in our DMF theory [Eq. (20)-(27)] (N eq = 2NN tr in direct simulations where N tr denotes 
the number of trials). In the case of iV = 100, for example, we get A* r e<? =20300 in the RT 
method, which is much larger than N eq = 8 in DMF theory. Our DMF theory successfully 
reduces the number of DEs, by taking account p K , 7 Kj a and p Kj A for global variables as well 
as local variables instead of m l K and for local variables. Although our DMF theory 
neglects spatial fluctuations in state variables, it has advantages of a tractable small number 
of DEs and clear semi-analytical nature, from which some qualitative results may be deduced 
without numerical calculations, as will be shown shortly [e.g. Eq. (52)]. When couplings 
have the spatial dependence as w — > Wij in Eq. (3), we have to rely on Eqs. (37)-(41) in the 
moment method. 



C. Property of DMF theory 

In this subsection, we will discuss the property of our DMF theory. It is possible to 
regard DMF theory as the single-site mean-field theory. Let us assume a configuration in 
which a single neuron % is embedded in an effective medium whose effect is realized by a 
given neuron i as its effective external input through the coupling w. We replace quantities 
of m*, C*'* and (1/iV) Z)fc(^) C«'a 111 coupling terms of Eqs. (37)-(41) by effective quantities 
of fM K , 7 Kj a and p Ki \ — (l/N)y Kt \, respectively. Then, in order to determine these quantities 
just introduced, we impose the single-site mean-field conditions given by [see Eqs. (46)-(48)] 

(49) 
(50) 

I E C& (51) 

Note that Eqs. (49)-(51) are assumed to hold independently of i and that m\ and C l ^\ in their 
righthand sides are functions of p K , 7 Kj a and p K; \. Conditions given by Eqs. (49)-(51) yield 
the self-consistent DEs for 7 k a and p K ^\ which are again given by Eqs. (20)- (27). The 



Pk = 

1k,X = 

1 

Pk,X T77k,A — 
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single-site approximation given by Eq. (49)-(51), which assumes that the quantities averaged 
at a given site are the same as those of the effective medium, is common in mean-field theories 
such as the Weiss theory for magnetism WJ§ and the coherent-potential approximation for 
random alloys [[5^ . 

We should note that the noise contribution is j3 2 in Eq. (22) while that is j3 2 /N in Eq. 
(25). It is easy to see that in the case of no couplings, we get 

P«,a = ^jf, (for w = 0) (52) 

which agrees with the central-limit theorem. On the other hand, in the case of (3 = 
and w 7^ 0, we get p K: \ = ^ K> \. Thus the ratio: p Kj a/7k,a changes as model parameters 
are changed. We will show that these changes in p K ^and A reflect on the firing time 
distribution and the degree of synchronous firings in neurons ensembles. 

Firing-Time Distribution 

The (nth) firing time of a given neuron % in the ensemble is defined as the time when 
Xi(t) crosses the threshold 9 from below: 

toin = {t | Xi(t) = 9; ±i > 0; t > t oin _ x + r r }, (53) 

where r r denotes the refractory period introduced so as to avoid multiple firings in a short 
period arising from fluctuations in voltage variables around the threshold. We get the 
distribution for the membrane potential variable x$ given by (for details see Appendix C) 

P(xi) * (-) <K^^), (54) 

where <f>(x) is the normal distribution function given by 

1 x 2 



with 

°i = VliA- (56) 

This implies that the distribution of the voltage variable Xi(t) is described by the Gaussian 
distribution with the mean of fii(t) and the variance of 71,1 (t). The probability given by Eq. 
(54) depends on the time because of the time dependence of Xi(t) and ag(t). The probability 
Woi(t) when Xi(t) at t is above the threshold 9 is given by E| 



W ol (t) = l~i;( 9 —^), (57) 
ere 

where ip{y) is the error function given by integrating <p(x) from —00 to y. Then the proba- 
bility averaged over the ensemble is given by 

i 

= 1 " <A(^)- (58) 

0-£ 
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The fraction of a given neuron i emitting output spikes at t is given by 

Z((t) = ^e Wl ) = ^)^)e te ). (59) 

at o i at oil 

where \i\ = d/ii/dt. When we expand /ii(t) in Eq. (59) around t* where H\{t*) = 9, it 
becomes 

W~*(^)4(£)e(,ii), (60) 

dt oe dt a e 



with 



Stot = —, (61) 



where /ii, /ii and are evaluated at t = i*. Zg(t) provides the distribution of firing times, 
showing that most of firing times of neurons locate in the range given as 

t eK- Stot, t* + 5t oe }. (62) 

In the limit of vanishing /?, Eq. (60) reduces to 

Z,(f) =*(*-*:)■ (63) 

Similarly, we get the distribution for the global variable X given by (for details see 
Appendix D) 

P(X) ~ (1) 0(^^), (64) 

with 

Og = VPlil- ( 65 ) 

This implies that the distribution of global voltage variable X(t) is described by the Gaussian 
distribution with the mean of fJ>i(t) and the variance of Pi,i(t). If we define the rath firing 
time relevant to the global variable X(t) as 

t gm = {t | X(t) = 9; X(t) >0;t> t gm ^ + r r }, (66) 

the fraction of firing around t — t* is given by 

Z g (t) = 0(^) ^) 9^), (67) 

0t O9 at (Jg 



with 



ot os = % (68) 
/ii 



Then most of t og locate in the range given by 
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t g e[t* -5t 09 , t* + 6t og }. (69) 

Since p^i is generally smaller than 7^1, we get a g < ae and 5t g < 5t g. In particular, in the 
case of no couplings, Eqs. (52), (56), (61), (65) and (68) lead to 

5t og = ^L. (for w = 0) (70) 

VN 



Synchronous Response 

Now we consider the quantity given by 

R(t) = ^2 E < M*) - X M 2 >= 2(7i,i - Pi,i). (71) 

When all neurons are in the completely synchronous state, we get Xi(t) = X(t) for all i, 
and then R(t) = 0. On the contrary, in the asynchronous (random) state, we get R(t) = 
2(1 — 1/N)~f ltl = Ro(t). Then the quantity defined by 

s{t) = 1 - myMt) = ip %~'^( N \ («) 

is 1 for the completely synchronous state and for the asynchronous state. We hereafter 
call S(t) the synchronization ratio, which provides the degree of synchronous firings in the 
ensemble. We get S = for (3 7^ with no couplings (w — 0), and S = 1 for w 7^ with no 
noises (f3 — 0). 



D. Calculated results 

We expect that our DMF equations given by Eqs. (20)- (27) may show bifurcation, 
synchronous and asynchronous states as well as chaotic states. In this study, we pay our 
attention to the response of the FN neuron ensembles to a single spike input, I^(t) given by 
Eq. (4), which is applied to all neurons in the ensemble. We have adopted the parameters 
of 9 = 0.5, a = 0.1, r r = 10, A = 0.10, t in = 100 and T w = 10. Parameter values of w, (3 
and N will be explained shortly. We get the critical magnitude of A c = 0.0442 below which 
firings of neuron defined by Eq. (53) cannot take place without noises {(3 = 0). We have 
adopted the value of A = 0.10 (> A c ) for a study of the response to a supra-threshold input, 
related discussion being given in Sec. IV. 

DMF calculations have been made by solving Eqs. (20)-(27) by the forth-order Runge- 
Kutta method with a time step of 0.01. Direct simulations have been performed by solving 
27V" differential equations as given by Eqs. (1) and (2) by using also the forth-order Runge- 
Kutta method with a time step of 0.01. Simulation results are the average of 100 trials 
otherwise noticed. Initial values of variables are set to be Hi = \xi = 7^1 = 7^1 = 7^2 = 
Pi,i — P2,2 = Pi, 2 = in DMF calculations, and xi = = for i = 1 to N in direct 
simulations. All calculated quantities are dimensionless. 

The time courses of means of \i\ and /i 2 calculated with f3 = 0.01, w = 0.0 and N = 100 
are shown in Figs. 1(a) and 1(b), respectively, where solid curves denote the results of DMF 
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theory and dashed curves those of direct simulations. We note that [i\ and [i^ obtained by 
two methods are in very good agreement and they are indistinguishable. At the bottom of 
Fig. 1(a) an input spike is plotted [see also Fig.2(a)]. States of neurons in an ensemble when 
an input spike is injected at t = 100, are randomized because noises have been already added 
since t = 0. Figures l(c)-l(h) show the time courses of various variances and covariances. 
Agreements between the two methods are good for 7^1, p^i, 7^2 and p\^- There is a fairly 
good agreement for 72,2 and p2,2- Comparing Figs. 1(c), 1(e) and 1(g) to 1(d), 1(f) and 1(h), 
respectively, we note that the relation given by Eq. (52): p K A = 7 Kj a/100 valid for w = 0, is 
supported by simulations: note that results in Figs. 1(d), 1(f) and 1(h) are multiplied by a 
factor of hundred. 

Figures 2(a) shows a single spike input, which is applied at t = 100 with a duration 
of T w = 10. The solid curve in Fig. 2(b) express Zg, the firing probability of the local 
variable Xi(t), which is a positive derivative of Wi shown by the dashed curve [Eqs. (58) 
and (59)]. They are calculated for f3 = 0.01, w = 0.0 and N = 100 in DMF theory. For a 
comparison, the simulation result for Z e is plotted in Fig. 2(c). Firings of neurons occur at 
t Q ~ 104 to 105 with a delay of about 4 ~ 5. Fluctuations of firing times of local variables 
5t e are 0.37 calculated by Eq. (61) in DMF theory, and 0.41 in simulations which is the 
root-mean-square (RMS) value of firing times defined by Eq. (53). In contrast, dashed and 
solid curves in Fig. 2(d) show W g and Z g , respectively, for the global variable X(t) in DMF 
theory, while Fig. 2(e) shows Z g obtained in simulations. Fluctuations in spike timings of 
the global variable are 5t og = 0.037 calculated by Eq. (68) in DMF theory, and 0.041 in 
simulations which is the RMS value of firing times defined by Eq. (66). These figures of 
5t og for the global variable are ten times smaller than respective values of 8t £ for the local 
variable. 

Noise-strength {(3) dependence 

We expect that as the noise strength is more increased, the distribution of membrane 
potentials is more widen and fluctuations of firing times are more increased. Filled squares in 
Fig. 3(a) show the (3 dependence of St Q e obtained by DMF theory [Eq. (61)] with w = 0.0 and 
N = 100, while open squares express the RMS value of firing times obtained by simulations. 
The agreement between the two methods is fairly good. In contrast, filled circles in Fig. 
3(a) show the (3 dependence of St og relevant to the global variable obtained by DMF theory 
[Eq. (68)] and open circles stand for RMS values of firing times in simulations. We note 
that 5t og is much smaller than 5t e although both 5t og and St e are proportional to (3 for 
weak noises under consideration. Figure 3(b) will be explained shortly in connection to the 
result of the w dependence. 

Cluster-size (N) dependence 

Filled squares in Fig. 4(a) show the N dependence of the local fluctuation of 5t Q e for 
(3 = 0.01 and w = 0.0, obtained by DMF theory, while open squares express that obtained 
by simulations. Simulations have not been performed for N > 100 because of a limitation 
in our computer facility. We note that 5t Q £ is independent of N because of no couplings 
(w = 0). In contrast, filled circles in Fig. 4(a) show the N dependence of the global 
fluctuation of 5t og obtained by DMF theory while open circles that by simulations. The 
relation: 5t og oc (l/y/N), holds as given by Eq. (70). Figure 4(b) for finite w will be 
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discussed shortly. 

Coupling-strength (w) dependence 

So far we have neglected coupling w among neurons, which is now introduced. Filled 
squares in Fig. 3(b) show the (3 dependence of local fluctuations of bt ol calculated by DMF 
theory for w = 0.2 and N = 100, while open squares that obtained by simulations. Filled 
and open circles express global fluctuations of 5t og in the DMF theory and simulations, 
respectively. Comparing these results with those for w = 0.0 shown in Fig. 3(a), we note 
that 5t Q £ is much reduced as w is increased although there is little change in St og . 

This is more clearly seen in Fig. 5(a), which shows the w dependence of firing-time 
fluctuations. Filled squares in Fig. 5(a) show fluctuations of 5t Q £ for the local variable 
obtained for (3 = 0.01 and iV = 100 by the DMF theory while open squares express those 
calculated by simulations. Filled and open circles in Fig. 5(a) show fluctuations of 5t og for 
the global variable obtained by the DMF theory and simulations, respectively, When w is 
increased, 5t Q £ is considerably decreased whereas 5t og is almost constant. Figure 5(b) shows 
a similar plot of the w dependence of firing times when the size of an ensemble is reduced 
to iV = 10. We note that 5t og for iV = 10 is 3.16 times larger than St og for N = 100 because 
5t og is proportional to 1/ \fN. 

Results obtained by DMF theory are analyzed in the Appendix E, where we get the 
expression for the w- and iV-dependent 5t Q £ given by [see Eq. (El)] 



where n — 1, 5t O £(0, 1) = 2.71, a x = 7.0 and a2 = —11.0. Bold, dashed curves for w < 0.2 
in Figs. 5(a) and 5(b) show the w dependence of 5t a p for N = 100 and 10, respectively, 
expressed by Eq. (73), which are in good agreement with results of DMF theory shown by 
filled squares. 

Log-log plots of Fig. 4(b) show the N dependence of 8t £ (squares) and 5t og (circles) for 
w = 0.2 and iV = 100, filled and open symbols denoting results of DMF and simulations, 
respectively. Although 5t og oc (l/y/N) as in the case of w = [Fig. 4(a)], 5t Q £ shows the 
peculiar N dependence, which arises from the (1 — 1/N) term in Eq. (73). The N dependence 
given by Eq. (73) with n — 1 and 2 is shown by thin solid curves at the uppermost in Fig. 
4(b). The result with n = 1 is in better agreement with the result of DMF theory shown by 
small filled squares than that with n = 2 (see Appendix E). 

Couplings among neurons work to increase the synchronous dynamics and to suppress 
local fluctuations. Figures 6(a) and 6(b) show the time sequence of the synchronization 
ratio S(t) defined by Eq. (72) for w — 0.1 and w = 0.2, respectively, with (3 = 0.01 and 
N = 100. Solid and dashed curves in Figs. 6(a) and 6(b) show results in DMF theory and 
simulations, respectively. Both results are in fairly good agreement. We note that S(t) has 
two peaks at times when pi,i(i) also has double peaks [Fig. 1(d)]. The maximum value of 
S(t) for w = 0.2 is S' max =0.132, which is larger than S max — 0.041 for w = 0.1. This trend 
is more clearly seen in Fig. 7 where the maximum magnitude of S(t), S max , is plotted as a 
function of w for iV = 10, 20, 50 and 100. It is shown that S max is increased as the coupling 
strength is increased, as expected. Figure 7 also shows that the effect of coupling is more 
significant in ensembles with smaller N. 




(73) 
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An analysis of the result obtained in DMF theory yields the expression for w- and N- 
dependent S max given by [see Eq. (E8) in Appendix E] 



S ma x(w, N) =CiW + C 2 W + .., (74) 



with 



«. = 4)<1 - 1) K (75) 

C2 = (i)(l-i)[62 + (l-i) 2 6?], (76) 

where b\ = 22 and b 2 = —290. Bold, dashed curves for w < 0.2 in Fig. 7 show S max 
expressed by Eqs. (74)-(76), which are in good agreement with results obtained in DMF 
theory shown by solid curves. If we define the coupling constant w m (N) for which S max is, 
for example, 0.3 for a given N, we get w m (A)=0.101, 0.147 0.237 and 0.322 for N=10, 20, 50 
and 100, respectively, which lead to w m (N)/w m (10)= 1.0, 1.46, 2.35 and 3.19, respectively, 
when w m (N) is normalized by w m (N = 10). This suggests that we may get w m (N) oc y/N. 
This arises from the fact that the relation: S max oc w 2 /N is nearly hold for S max = 0.3 for 
which the contribution from the w 2 term is more considerable than that from the w term 
in Eqs. (74)- (76). Of course, it is not the case for much smaller value of S max for which the 
first term is more dominant than the second term. 

Expressions of Eqs. (E1)-(E8) for w- and A-dependence of fluctuations and the synchro- 
nization ratio, which are obtained based on the results calculated in DMF theory, are useful 
in a phenomenological sense. For example, in the case of negative (inhibitory) couplings, 
Eqs. (El) and (E8) yield an increase in St e and a negative S, which are supported by numer- 
ical calculations with DMF theory and simulations (not shown). We have tried to extract 
coefficients ai, a 2 , b\ and b 2 in Eqs. (E1)-(E8), by expanding Eqs. (20)-(29) in terms of w, 
but have not succeeded yet. 



III. LARGE CLUSTER CONSISTING OF MULTIPLE SUB-CLUSTERS 

A. Formulation 

It is possible to extend our DMF theory to a large FN neuron cluster which is divided 
into multiple M sub-clusters according to their functions. Dynamics of a single FN neuron 
j in a given sub-cluster m (=1 to M) which consists of N m neurons, is described by the 
nonlinear DEs given by 

F[ Xi (t)) ~ cy t {t) + 4 cl) (t) + 4 c2) (t) + lt\t) + CM, (77) 
bxi(t) - d Vi (jt) + e, (i = l-N m ) (78) 

where Xj and yi denote the fast (voltage) variable and slow (recovery) variable, respectively, 
the Gaussian white noise with < £j(t) >= and < £j(t) £/(£') >= 01 8(t — t'), the 
bracket < ■ > denoting the average In Eq. (77), l\ if) and ij; (t) given by 



dxi(t) 

dt 
dyjjt) 
dt 
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4 cl) (t) = (w mm /N m ) E G(xj(t)), (79) 

j em 

lf\t)= E («W^») J2G(x k (t)), (80) 

express the couplings within the sub-cluster m with the strength w mm , and those between 
sub-clusters with the strength w mn , respectively, N rn the number of neurons in the sub- 
cluster m, and G(x) is the sigmoid function. I$(t) stands for an external single spike input 
applied to all neurons in the sub-cluster m, as given by Eq. (4). 

As in the Sec. IIA, we first define the global variables for the sub-cluster m by 



X m (t) = (1/N m ) £ xi(t), (81) 
Y m (t) = (1/N m ) E Vi(t), (82) 



and their averages by 

A*?*(*) = < X rn {t) >, (83) 
*V(t) = <Y m (t)>, (84) 

Next we define variances and covariances between intra-neuron variables, given by 

iZ = E < ( 6x T) 2 >, (85) 

7 2 m 2 = ^ E < W) 2 >, (86) 
7i? 2 = ^ E < W V) >, (87) 

where 5x™ = Xi(t)—pL m (t) and 5y™ = Di{t)— n™(t), and those between inter- neuron variables, 
given by 

p™ = < (SX m ) 2 >, (88) 

p™ 2 = < (5Y m ) 2 >, (89) 
p™ = < (5X m 5Y m ) >, (90) 

where 5X rn = X m {t) - pi m (t) and 5Y m = Y m {t) - (i™(t). After some manipulations, we get 
the following differential equations: 

du m 1 

= /o" 1 + ^l?! - C/# + «W (1 - — ) U om + E Wmn ^On + /<?(*), (91) 
flr iVm n(^m) 

^ = b^T- d& + e, (92) 
= 2[(/r + 3/ 3 - 7 -) 7 ri - c 7l m 2 ] + 2w mm (p? A - j^) C/ lm 
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+2 Yl WmnPltUm + Pl, (93) 

di m 

-j^ = 2(b 1 ™-d 1 ™), (94) 
d^ m T m 

—jj- = 6 7l,l + Ul + 3/ 3 7l,l - «)7l,2 - C 7 2 ,2 + ^mm(Pl,2 ~ J^) U\m 

U ln , (95) 

do m 1 /? 2 

-gi = 2[(/r + 3/ 3 m 7 - )p™ - cp™ ] + 2 Wmm (1 - — ) p™ tf lm + ^, (96) 

J^ = 2(bp? !2 -dp™ 2 ), (97) 

do m 1 

= + if? + 3/3 m 7i m i - d)pl 2 - cp% 2 + w mm (1 - — ) p™ U lm , (98) 

with 

f/ 0m = C + ^7ri, (99) 

u lm = gT + 3g?iZ> ( 10 °) 

where /f = (l/£\)F^(jif), = (1/£\)G®(ji?). and (3^ = ^ Ei Gm A 2 - Now we have to 
solve 8M dimensional deterministic differential equations, which is more amenable than to 
solve 2NM stochastic DEs. 

When a given cluster can be divided into excitatory and inhibitory sub-groups and when 
variance and covariance terms in Eqs. (91)-(98) are neglected, we get 

— -fo- c/if + w EE U E + w EI U Z + 4 e) (t), (101) 



dt 

E 

2 -bpf-dp e 2 + e, (102) 



dp E 



dt 
dp[ 
~dt 
dp{ 

dt 



f Q - cp{ + w n Uj + w IE U E + Ij e \t), (103) 
bp[ - dp 2 + e, (104) 



where the supra-script E and / stand for the excitatory and inhibitory clusters, respectively. 
This corresponds to the result of Wilson and Cowan ||60|| . Then our DMF theory given 
by Eqs. (91)-(98) may be regarded as a generalized version of the Wilson-Cowan theory 
including fluctuations of local and global variables. 

Firing-Time Distribution 

The fraction of firings of neurons in the sub-cluster m is given by 

z om {t) = 0(^^) -(^M e(/ir), (los) 

(j£ m at ai m 

with 
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aim = yf^i- (106) 

When we expand fi m (t) in Eq. (105) around t* om where n m (t* m ) = 9, it becomes 

t — t* rl n m 

Z om (t) ~ <P(-^) -ri—) ©(/O, (107) 
01 om ox o~e m 



with 



Km = (108) 



mi 



where fi m , \i\ m and ai m are evaluated at t — t* om . This shows that most of firing times of a 
given sub-cluster m locate in the range given as 

torn £ \Pom St om , t om -\~ 5t orn j. (109) 



Synchronous Response 

The synchronization ratio of a given sub-cluster m is given by 

M) ~ (i-i/iv w ) ' (110) 

which is and 1 for completely asynchronous and synchronous states, respectively. 



B. Calculated results 

We have performed model calculations, assuming M (= 10) sub-clusters, each of which 
consists of N m (= N) neurons. They are connected by feed-forward inter-sub-cluster cou- 
pling given by w mn = w 2 S nm -x, which is allowed to be different from the intra-sub-cluster 
coupling given by w mm = wi for all m. A single spike input given by Eq. (4) is applied only 
to the first sub-cluster (m = 1), and an output of a sub-cluster m is subsequently forwarded 
to the next sub-cluster m + 1. This is conceptually similar to the synfire chain f6~i~| . When 



W2 is too small, signals cannot propagate through sub-clusters. The critical value of the 
inter-sub-cluster coupling w 2c , below which a spike cannot propagate through sub-clusters, 
is w 2c =0.064, 0.028 and 0.020 for io x =0.0, 0.1 and 0.2, respectively, with (3 = 0.0 and 
N = 100. 

Figure 8(a) shows the time course of Z om (t) calculated in DMF theory with f3 = 0.05, 
wi = W2 = 0.1 and iV = 100. Signals propagate through sub-clusters with 5t om ~ 0.8 
for all m. The result is in good agreement with that obtained in direct simulations (not 
shown). Synchronization ratios S m (t) shown in Fig. 8(b) have double peaks [see Figs. 6(a) 
and 6(b)]. The maximum value of S m (t), for example, is 0.022 for m = 1 at t — 122.2. 

In contrast, Fig. 8(c) shows the time course of Z om (t) for the increased noise intensity of 
(3 = 0.23, which shows that signals cannot propagate, dying out at the sixth sub-cluster. In 
this case, the agreement of DMF results with simulations is not satisfactory. Synchronization 
ratios S m (t) for (3 = 0.23 shown in Fig. 8(d) have multiple peaks for 1 < m < 4, double 
peaks for m = 5, a single peak for m = 6, and it disappears for m > 6. 
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Figure 9(a) shows the m-dependence of local fluctuations St om for various (3 with w\ = 
W2 = 0.1 and N = 100. We note that 5t om is almost constant for (3=0.05 and 0.10. In the 
case of (3 = 0.23, however, d~t om is divergently increased at m = 5. This behavior is not 
changed when we adopt a different set of parameters. Fig. 9(b) shows a similar plot of St om 
as a function of m for w% = 0.0, W2 = 0.1 and iV = 100. Signals propagate with 5t om = 
0.04 and 0.12 for [3=0.01 and 0.05, respectively. For (3 = 0.09, however, a spike dies out at 
m = 8. 

Figure 10 shows the W\ dependence of the critical noise strength f3 c above which signals 
cannot propagate. We get (3 C = 0.09 and 0.23 for W\ = 0.0 and 0.1, respectively, for N = 100 
as discussed above. When w\ is set to be 0.2, (3 C becomes 0.38 for iV = 100. We note that 
(3 C is almost linearly increased by increasing w±. Figure 10 also shows that the critical value 
of (3 C becomes larger as the size of sub-cluster (N) is larger. 

IV. CONCLUSION AND DISCUSSION 

We have proposed DMF theory for stochastic FN neuron ensembles, in which means, 
variances and covariances of local and global variables are taken into account. DMF theory 
has been shown to be derived in various ways: series expansions of means, variances and 
covariances of local and global variables (Sec. IIA), re- arrangement of moments in RT's 
method (Sec. IIB) and a single-site approximation to RT's method [Sec. IIC] [Q. Our 
DMF theory, which assumes weak noises and the Gaussian distribution of state variables, 
goes beyond the weak coupling because no constraints are imposed on the coupling strength. 
Calculated results based on DMF theory are in fairly good agreement with those obtained 
by direct simulations for weak noises. When the noise intensity becomes stronger, the state- 
variable distribution more deviates from the Gaussian form (see Fig. 3 of |47[j), and the 
agreement of results of DMF theory with those of simulations becomes worse. Nevertheless, 
our DMF theory is expected to be meaningful for qualitative or semi-quantitative discus- 
sion on the properties of neuron ensembles or clusters. It is possible to regard nonlinear 
differential equations given by Eqs. (20)-(27) [or Eqs. (91)-(98)] as the mean-field FN model 
for neuron ensembles or clusters. We hope that our DMF theory may play a role of the 
molecular- field (Weiss) theory in magnetism ||57|| : the Weiss theory provides a clear phys- 
ical picture on various magnetic properties despite some disadvantages such that it yields 
too-high critical (Curie) temperature, wrong critical indices and wrong temperature depen- 
dence for magnetization at low temperatures. Our DMF theory may be applied to a general 
conductance-based nonlinear systems. When it is applied to an ensemble of iV HH neurons, 
we get the 24 deterministic nonlinear differential equations, which are more amenable than 
original 4iV stochastic equations. Furthermore our DMF theory based on means, variances 
and covarinces of local and global variables can be applied to more general stochastic systems 
besides neural networks. 

In summary, we have developed a semi-analytical DMF theory for FN neuron ensembles. 
In order to show the feasibility of the DMF theory, we have studied the response of ensembles 
of FN neurons to a single spike input. The result is summarized as follows: (i) the spike 
timing precision of the global variable is much improved by increasing the ensemble size, 
even when there is no couplings among constituent neurons, (ii) by increasing the coupling 
strength, the spike transmission is enhanced by the synchronous response, and (iii) the 
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spike propagation with a fairly precise timing is possible in large-scale clusters when the 
noise strength is moderate. The origin of the item (i) is the same as that yielding the 
central-limit theorem. Couplings work to suppress local fluctuations and to increase the 
synchronization ratio [Eq. (72)]. Items (i) and (ii) are consistent with the results reported 
previously ||30|| - f35|j . The item (iii) agrees with the result of recent simulations for synfire 
chains, each layer of which consists of 100 IF neurons |62j| . Items (i)-(iii) are beneficial to the 
population temporal-code hypothesis mentioned in the introduction. Although calculations 
reported in this paper have been limited to supra-threshold inputs, it is possible to study 
the response to sub-threshold inputs with the use of DMF theory. We may investigate 
combined effects of white noises and the heterogeneity in model parameters, which have 
been intensively studied in recent years [ S3| . Such calculations are in progress and will be 
reported in a separate paper. 
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APPENDIX A: DERIVATION OF EQS. (20)-(27) AND EQS. (32)-(36) 

From Eqs. (9)-(12), we get the differential equations for the deviations of 5xi and 5yi of 
the neuron i, given by 

d° Xi " - ■ » "• „2 _ \ j f s„3 _ „ X „, , t J AfW 



dt 

dSyi 
dt 

with 



hSxi + Wx\ - 71,1) + fM - c5 Vi + & + 5I?>, (Al) 
bSxi — d5yi, (A2) 



*4 C) = w | E Sxj + gA E - (1 - ^)7i,x] + | E ^ I (A3) 
The differential equations for the variances and covariances are given by 

-^f = ^^E < ^i^i + (, Sx t s Vi) ^i^A2 + (<fy 2 ) 5 k2 5 X2 ] >, 

= i E < {2[^ (^*)] 5 Kl 5 xl + [S Vi (^) + 8* (^)l *«ifc> 
+ 2[6y t (^)] 6 k2 6 X2 } >, (A4) 

= ^(V^ 2 ) E E < [(SxiSxj) 6 Kl 5 xl + (5xi5yj) 5 Kl 5 X2 + (Sy^Vj) S k2 5 X2 ] > 

i 3 
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= ^ E E < {2[&* (^)l Wai + [<fo (^) + A* (^)]«a 2 

* 3 

+ 2[5 Vi (^)] 5 K2 5 A2 } >, (A5) 

In the process of the calculation using Eqs. (20)-(27), we have adopted the following ap- 
proximations: 



(1) the forth-order variances are assumed to be [47 



i < Sx i > = 3 7i,i7i,i, ( A6 ) 

i 

— < Sx 3 Syi > = 3 7i,i7i,2, (A7) 



and 

N 2 
1 

iV2 



4EE< > = 3 7i,iPi,i, (A8) 

3 

' E E < 6 Vi 6x ) > = 3 7i,iPi,2, (A9) 

* 3 



other forth terms being set zero. 

(2) the third-order variances and terms higher than forth order are neglected, 

Calculations of Eqs. (32)-(36) are siimilar to those discussed above if we read as m\ — > fix, 
m 2 ~^ ^2, Ax; —> and Ayi — > by^. 



APPENDIX B: SOME LIMITING CASES OF EQS. (20)-(27) 

(1) In the limit of a single (N = 1) neuron, Eqs. (20)-(27) reduce to 



^ = /o + / 2 7i,i-c/x 1 + / (e) (t), (Bl) 

— — = bui - dfi 2 + e, (B2) 

^ = 2(/ lTl)1 + 3/ 3 7i 2 1 - C7i, 2 ) + P 2 , (B3) 

^ = 2(6 7l , 2 -rf 7 2,2), (B4) 

'-TT = &7i,i + (/i - d)7i,2 + 3/ 3 7i,i7i,2 - C72, 2 , (B5) 

p K ,A = 7«,A- (B6) 



Equations (B1)-(B5) agree with the results of Rodriguez and Tuckwell (RT) |44| and Tanabe 
and Pakdaman (TP) |47[]. In RT, the forth terms which appear in the process of calculating 
dji t i/dt and dji^/dt in Eqs. (B3) and (B5), are assumed to be zero, whereas in TP, they 
are assumed to be as given by Eqs. (A6) and (A7). 
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(2) In the limit of large N, where the exclusion of the self-couplings in Eq. (3) may be 
neglected, Eqs. (20)-(27) become 



dt 
~dT 

_ or. t o... ... i o , ? r> 

dt 

dl2,2 

dt 
dt 

^ = 2[(/i + 3/s 7 i,i)Pi,i - cpi, 2 ] + 2wp ltl U 1 + ^, (B12) 



fo + f 211,1 -cfi 2 + wU + I {e \t), (B7) 

bfxi - dfj,2 + e, (B8) 

2[(/i + 3/ 3 7i,i)7i,i - 071,2] + 2 WPlil C/i + (3\ (B9) 

2(671,2-^72,2), (BIO) 

67l,l + (/l + 3/371,1 - ^)7l,2 - C72,2 + wp lj2 U!, (Bll) 



2(6/01,2-^,2), (B13) 

6^1,1 + (/l + 3/371,1 - rf)pi j2 - Cp 2 ,2 + W/01,2 ^1, (B14) 



where U = g + 5271,1 and ^1 = #1 + 3^ 3 7i,i- 

APPENDIX C: DERIVATION OF EQ. (54) 

The distribution P{xj) in Eq. (54) is formally given by 

P(xi) = J...J n j tf i) dx j n j dy j p(x 1 ,....,x N ,y 1 ,....,y N ), (CI) 

with the probability distribution function (pdf) of p(xi, xn, Vi, Vn) f° r the 2N- 
dimensionl vector z = (x±, ....,Xn,Ui, ■■■■,Un), given by 

p(x 1 ,....,x N ,y 1 ,....,y N ) = 1 exp[~(z - nf V' 1 (z - p)], (C2) 

(2n)»J\V\ 2 

where /x and V express the mean vector and the variance-covariance matrix, respectively, 
In the case of a single FN neuron (N — 1), pdf is given by 

1 1 

p(xi,yi) =Pi(x 1 ,y 1 ) = j== exp[--(z - /*)* V" 1 (z -/*)], (C3) 

(27T) V VVn 2 

with 

z=(x 1 ,?/ 1 ) t , (C4) 
A* = Afc)*, (C5) 

V = f ^ 7l ' 2 ) . (C6) 

V 7l,2 72,2 / 
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Substituting Eqs. (C3)-(C6) to Eq. (CI), we get 



P(x 1 ) = / dyi p(x u yi) 




(C7) 



where 0(x) denotes the normal distribution function: 

4>{x) = -l=exp(-^). (C8) 

In the case of arbitrary TV under consideration, the calculation of P(xi) may be per- 
formed within DMF as folllows. As mentioned in Sec. IIC, our DMF theory assumes the 
configuration in which a single neuron is embedded in an effective medium characterized by 

7k,a and p Kt \ [Eqs. (49)-(51)]. Thus it is effectively the problem of a single neuron in the 
effective medium. Means (/i K ), variances (7 k ,a) and covariances (p K ,\) of local variables are 
determined by Eqs. (20)-(24). Then the calculation of P(xi) for N > 1 is the same as that 
for N = 1 mentioned above, and it is given by 



Equations (20), (21), (25)-(27) form DEs for means (fi K ), variances (7 K) a) and covariances 
(Pk,x) f° r global variables, X and Y. Then P(X) in Eq. (64) is given by 




(C9) 



APPENDIX D: DERIVATION OF EQ. (64) 




(Dl) 



with pdf for the two-dimensional vector z = (X, Y)* given by 




(D2) 



with 



fi = (//i,// 2 )*, 



(D3) 




(D4) 



Substituting Eqs. (D2)-(D4) to Eq. (Dl), we obtain 




(D5) 



where <f>(x) denotes the normal distribution [Eq. (C8)]. 
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Alternatively P(X) is expressed by 

P(X) = J ...J n i dx i n i dy i p(x 1 ,...,x N ,y 1 ,...,y N )5(X--jyYl x i)i ( D6 ) 

i 

where p(x±, ...,xn, yi, Vn) stands for pdf for 2iV-dimensional vector [Eq. (C2)]. However, 
a calculation of P(X) based on Eq. (D6) is difficult except for the no coupling case (w — 0), 
for which pdf is given by 

p(x u ...,x N , yi, ...,y N ) = RjPi(xi, yi), (D7) 

Pi(xi) being pdf for N — 1 [Eq. (C3)]. Performing integrals with respect to yi in Eq. (D6) 
with Eq. (D7), we get 



p(x) =/.../ n,-^ 0(^—^1 5(x-iE4 



(D8) 



By using the procedure conventionally used for proofing the central-limit theorem, we obtain 
Eq. (D5) with = r )i t \/\fN (for w = 0). We should note that, a calculation of P(X) 
based on Eq. (Dl) is easier than that based on Eq. (D6) and that the former is applicable 
for finite couplings. 



APPENDIX E: ANALYSIS OF NOISE, COUPLING AND SIZE DEPENDENCE 

(1) Sta and 5t og 

Based on the calculated results of DMF theory, we have tried to obtain the analytical 
expression of (3-, w- and iV-dependence of St a e and 8t og . Figures 3(a) and 3(b) show that 5t Q £ 
and St og are proportional to j3 for weak noises, for which both 7^1 and pi t i are proportional 
to (3 2 [see Eqs. (56), (61), (65) and (68)]. From results shown in Figs. 4 and 5, we have 
obtained expressions given by 

6UW ' N) - l-(h(l-lr(a lW + a^ + ..), (El) 



5t oe (0,l) 2 yv V 

5tog(W,N) 1 

5t oe (0, 1) ~ ViV' 



(E2) 



where n = 1, 5t o e(0, 1)=2.71, ai = 7.0 and a 2 = —11.0. The N dependence of 5t Q e expressed 
by Eq. (El) with n = 1 and 2 are shown by thin solid curves at the uppermost in Fig. 
4(b) with DMF result (small filled squares): these results are shifted upward by 0.433 for 
a clarity of the figure. The result with n = 1 is in better agreement with the DMF result 
than that with n = 2. On the other hand, bold, dashed curves in Fig. 5(a) and 5(b) show 
the w dependence of 5t e for iV = 100 and 10, respectively, expressed by Eq. (El) with 
n = 1, which is in good agreement with results of DMF theory shown by filled squares. This 
implies from Eqs. (56), (61), (65) and (68) that the w- and iV-dependence of 7^1 and p± t i 
evaluated at t — t* where ^i{t*) = 0, are given by 
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7i,i(™>W) i . 1. ^ , 



7i,i (0,1) 1 V 

Pi.iKAO 1 



1-(1-— ) ( ai w + a 2 w 2 + ..), (E3) 

(E4) 



7i,i(0,l) AT 

Note that 5t o e(0, 1) and 7i,i(0, 1) are proportional to (3 and /3 2 , respectively. 

(2) S max 

In order to discuss the expression of /?-, u>- and A^-dependent S^ax, we have analyzed 
results of S 1 ,™^ shown in Fig. 7 by 

Smaz = CiW + C 2 W 2 + (E5) 

to guess how expansion coefficients of c\ and c 2 depend on N. After several tries, we have 
concluded that the w- and A^-dependence of 71,1 and evaluated at t — where pi,i(t) 
has the maximum value, may be given by 

+ + (E6) 

7i,i(U, 1) A 
Pi,i(w,N) l_ 
7i,i(0, 1) ~iV' 

yielding given by [see Eq. (72)] 



(E7) 



S max (w, N) = (1)(1 - I)- 1 + [fe 2 + (1 - 1) 2 6 2 ] t, 2 }, (E8) 

where m = 2, &! = 22 and 6 2 = —290. Bold, dashed curves in Fig. 7 show the w dependence 
of S m ax expressed by Eq. (E8) for various A^ values, which are in fairly good agreement 
with results of DMF theory shown by solid curves. We should point out that a factor of 
(1 — l/N) in Eqs. (El), (E3), (E6) and (E8) appears because the coupling w does not work 
in a single-neuron case (N = 1) and that at least the second power (m = 2) is necessary 
in Eq. (E6) for S max to vanish in the N = 1 limit. A functional form of Eq. (E3) may be 
different from that of Eq. (E6) because the former is evaluated at t* while the latter at t^. 
Our DMF calculation shows that when (3 is increased for a fixed (finite) w value, S max is 
gradually decreased, although Eq. (E8) has no (5 dependence. This is due to contributions 
of 0(/3 4 ) to 7! ! and p ltl , which have been not included in the above discussion. 
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FIGURES 

FIG. 1. Time courses of means, variances and covariances calculated by DMF theory (solid 
curves) and simulations (dashed curves): (a) fi±, (b) fj,2, (c) 71,1, (d) p^i, (e) 72,2, (f) P2,2, (g) 
7i j2 and (h) p lj2 , for A = 0.10, f3 = 0.01, w = 0.0 and N = 100. Results of (d), (f) and (h) are 
multiplied by a factor of hundred. The chain curve at the bottom of (a) expresses a single input 
spike, 1^ in Eq. (4) [see also Fig. 2(a)]. 

FIG. 2. Time courses of (a) I^ e \ (b) Wn (the dashed curve) and Z^ (the solid curve) in DMF 
theory, (c) Z(i in simulations, (d) W g (the dashed curve) and Z g (the solid curve) in DMF theory, 
and (e) Z g in simulations, for A = 0.10, (3 = 0.01, w = 0.0 and N = 100. 

FIG. 3. (a) The (3 dependence of 5t Q e (squares), and 5t og (circles) for w = 0.0 and (b) that for 
w = 0.2 with N = 100, filled symbols denoting results in DMF theory and open symbols those in 
simulations. 

FIG. 4. Log-log plots of 5t £ (squares) and 5t og (circles) against N for (a) w = 0.0 and (b) 
w = 0.2, filled symbols denoting results in DMF theory and open symbols those in simulations. 
Shown at the uppermost in (b) are the DMF result (small, filled squares) and results with Eq. 
(El) with n = 1 and 2 (thin solid curves): they are shifted upward by 0.433 for a clarity of the 
figure (see text). 

FIG. 5. The w dependence of 5t Q £ (squares) and 5t og (circles) for (a) N = 100 and (b) iV = 10 
with [3 = 0.01, filled symbols denoting results in DMF theory and open symbols those in simula- 
tions. Bold, dashed curves for w < 0.2 express Eq. (73) (see text). 

FIG. 6. The time course of synchronization ratio S for (a) w = 0.1 and (b) w = 0.2 with 
[3 = 0.01 and N = 100, solid curve denoting results of DMF theory and dashed curve those of 
simulations. 

FIG. 7. The w dependence of the maximum of S, S max , for N = 10 (squares), N = 20 
(triangles), N = 50 (inverted triangles) and N = 100 (circles) with (3 = 0.01, filled symbols 
denoting results of DMF theory and open symbols those of simulation. Bold, dashed curves for 
w < 0.2 express Eq. (74) (see text). 

FIG. 8. Time courses of Z om for (a) (3 = 0.05 and (b)(3 = 0.23 with t«i = u>2 = 0.1, and time 
courses of S m for (c) (3 = 0.05 and (d)(3 = 0.23 with w\ = W2 = 0.1, calculated for N = 100 and 
M = 10 by DMF theory. 

FIG. 9. 5t om as a function of m for (a) w\ = W2 = 0.1 and (b) w\ = 0.1 and W2 = 0.0 with 
N = 100 and M = 10. 
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FIG. 10. The w\ dependence of (3 C for various N values with 
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